How to exploit data on Pangeo

How to exploit data on Pangeo#

Snow and Cloud Cover Example#

The executable content of this notebook is adapted from EO-College/cubes-and-clouds

Authors & Contributors#

Authors#

Michele Claus, Eurac Research (Italy), clausmichele (original author)
Pier Lorenzo Marasco, Provare LTD (UK), @pl-marasco (author of the conversion)

Contributors#

Alejandro Coca-Castro, The Alan Turing Institure, acocac (author)
Anne Fouilloux, University of Oslo (Norway), @annefou
Justus Magin, UMR-LOPS CNRS(France), @justusmagin
Tina Odaka, UMR-LOPS Ifremer (France), @tinaok

Introduction#

In this exercise, we will build a complete the same EO workflow as OpenEO using cloud provided data (STAC Catalogue), processing it locally; from data access to obtaining the result.

We are going to follow these steps in our analysis:

Load satellite collections
Specify the spatial, temporal extents and the features we are interested in
Process the satellite data to retrieve snow cover information
Aggregate information in data cubes
Visualize and analyse the results

Relevant resources#

More information on Pangeo can be found here: https://pangeo.io/
More information on the STAC specification can be found here: https://stacspec.org/

Import libraries#

# Data Manipulation and Analysis Libraries
import pandas as pd  
import numpy as np 

# Geospatial Data Handling Libraries
import geopandas as gpd 
from shapely.geometry import mapping  
import pyproj

# Multidimensional and Satellite Data Libraries
import xarray as xr 
import rioxarray as rio
import stackstac

# Data Visualization Libraries
import holoviews as hv
import hvplot.xarray
import hvplot.pandas

# Data parallelization and distributed computing libraries
import dask
from dask.distributed import Client, progress, LocalCluster

# STAC Catalogue Libraries
import pystac_client

Connect to the Dask cluster#

Here we creates a Dask client, which is essential for managing and executing parallel computations efficiently in the subsequent parts of the notebook. There are situation where you can connect to a Dask Gateway, but for this exercise we will use a local cluster.

cluster = LocalCluster()
client = Client(cluster)

Client address can be copy and pasted to the dashboard to monitor the progress of the computations.

client

Client

Client-d516f08c-7a6d-11ee-8f59-6045bd836204

Connection method: Cluster object	Cluster type: distributed.LocalCluster
Dashboard: http://127.0.0.1:8787/status

Cluster Info

LocalCluster

a4dbacc9

Dashboard: http://127.0.0.1:8787/status	Workers: 4
Total threads: 4	Total memory: 15.61 GiB
Status: running	Using processes: True

Scheduler Info

Scheduler

Scheduler-eeb7764b-a71b-4e87-98d8-99e410646f9f

Comm: tcp://127.0.0.1:45411	Workers: 4
Dashboard: http://127.0.0.1:8787/status	Total threads: 4
Started: Just now	Total memory: 15.61 GiB

Workers

Worker: 0

Comm: tcp://127.0.0.1:37705	Total threads: 1
Dashboard: http://127.0.0.1:42115/status	Memory: 3.90 GiB
Nanny: tcp://127.0.0.1:34287
Local directory: /tmp/dask-scratch-space/worker-bhnn81yw

Worker: 1

Comm: tcp://127.0.0.1:36073	Total threads: 1
Dashboard: http://127.0.0.1:43909/status	Memory: 3.90 GiB
Nanny: tcp://127.0.0.1:37781
Local directory: /tmp/dask-scratch-space/worker-2sdgbscd

Worker: 2

Comm: tcp://127.0.0.1:34715	Total threads: 1
Dashboard: http://127.0.0.1:45119/status	Memory: 3.90 GiB
Nanny: tcp://127.0.0.1:40853
Local directory: /tmp/dask-scratch-space/worker-sf1j_l70

Worker: 3

Comm: tcp://127.0.0.1:45561	Total threads: 1
Dashboard: http://127.0.0.1:41387/status	Memory: 3.90 GiB
Nanny: tcp://127.0.0.1:44965
Local directory: /tmp/dask-scratch-space/worker-d0bakp5d

Load data#

Catchment outline#

We will use the catchment as our area of interest (AOI) for the analysis. The catchment is defined by a polygon, which we will load from a GeoJSON file. The GeoJSON file contains the geometry of the catchment in the WGS84 coordinate reference system (EPSG:4326) and that has to be defined.

aoi = gpd.read_file('../data/catchment_outline.geojson', crs="EPGS:4326")
aoi_geojson = mapping(aoi.iloc[0].geometry)

Satellite collections#

Search for satellite data using STAC#

We will utilize the STAC API to search for satellite data in this exercise, specifically leveraging the API provided by AWS/Element84. The STAC API operates as a RESTful service, enabling the querying of satellite data with various filters such as spatial range, time period, and other specific metadata. This API is constructed based on the STAC specification, a collaborative, community-driven standard aimed at enhancing the discoverability and usability of satellite data. Numerous data providers, including AWS, Google Earth Engine, and Planet (Copernicus Data Space Ecosystem (CDSE) is coming soon**), among others, have implemented the STAC API, exemplifying its widespread adoption and utility in accessing diverse satellite datasets.

For the purposes of this exercise, we will limit the search to the Sentinel 2 L2A collection, which is a collection of Sentinel 2 data that has been processed to surface reflectance (Top Of Canopy). We will also limit the search to the time period between 1st February 2019 and 10th June 2019 and to the extent of the catchment.

** at the moment of writing the STAC catalog of the CDSE is not yet fully operational.

URL = "https://earth-search.aws.element84.com/v1"
catalog = pystac_client.Client.open(URL)
items = catalog.search(
    intersects=aoi_geojson,
    collections=["sentinel-2-l2a"],
    datetime="2019-02-01/2019-06-10"
).item_collection()
len(items)

Get bands information#

As the original data provides bands with different names than the original Sentinel 2 bands, we need to get the information about the bands.

# Get bands information
# selected_item = items[1]
# for key, asset in selected_item.assets.items():
#     print(f"{key}: {asset.title}")

Load data#

We will use the stackstac library to load the data. The stackstac library is a library that allows loading data from a STAC API into an xarray dataset. Here we will load the green and swir16 bands, which are the bands we will use to calculate the snow cover. We will also load the scl band, which is the scene classification layer, which we will use to mask out clouds. Spatial resolution of 20m is selected for the analysis. The data is loaded in chunks of 2048x2048 pixels.

Stackstac is not the only way to create a xarray dataset from a STAC API. Other libraries can be used, such as xpystac or odc.stac. The choice of the library depends on the use case and specific needs.

ds = stackstac.stack(items,
                    bounds_latlon=aoi.iloc[0].geometry.bounds,
                    resolution=20,
                    chunksize=2048,
                    assets=['green', 'swir16', 'scl'])

/home/runner/micromamba/envs/bids23/lib/python3.11/site-packages/stackstac/prepare.py:408: UserWarning: The argument 'infer_datetime_format' is deprecated and will be removed in a future version. A strict version of it is now the default, see https://pandas.pydata.org/pdeps/0004-consistent-to-datetime-parsing.html. You can safely remove this argument.
  times = pd.to_datetime(

Calculate snow cover#

We will calculate the Normalized Difference Snow Index (NDSI) to calculate the snow cover. The NDSI is calculated as the difference between the green and the swir16 bands divided by the sum of the green and the swir16 bands.

For a matter of clarity we will define the green and the swir16 bands as variables. Other approaches can be used to manage the data, but this is the one we will use in this exercise.

green = ds.sel(band='green')
swir = ds.sel(band='swir16')
scl = ds.sel(band='scl')

Let’s compute the NDSI and mask out the clouds.

ndsi = (green - swir) / (green + swir)

Dask allow to persist the data in memory, which is useful to speed up the computation. The persist method will load the data in memory and will keep it there until the end of the analysis. More explanation about persist, load and other dask related topics will be given in the next part of the exercise.

ndsi = ndsi.persist()

We will mask out the clouds, which are identified by the values 8 and 9 in the scene classification layer (scl). The scl contains information about the type of land cover. We will mask out the clouds, which are identified by the values 8 and 9 in the scl layer.

More detailed info can be found here: https://sentinel.esa.int/web/sentinel/technical-guides/sentinel-2-msi/level-2a/algorithm-overview

snow = xr.where((ndsi > 0.42) & ~np.isnan(ndsi), 1, ndsi)
snowmap = xr.where((snow <= 0.42) & ~np.isnan(snow), 0, snow)
# mask = (scl != 8) & (scl != 9) & (scl != 3) 
mask = np.logical_not(scl.isin([8, 9, 3]))  # more elegant but not sure about it from a teaching perspective
snow_cloud = xr.where(mask, snowmap, 2)

Process snow cover data#

Mask data#

As we are only interested to the snow cover in the catchment, we will mask out the data outside the catchment. To achieve it we need to convert the catchment geometry to the same coordinate reference system as the data. The data is in the UTM32N coordinate reference system (EPSG:32632).

aoi_utm32 = aoi.to_crs(epsg=32632)
geom_utm32 = aoi_utm32.iloc[0]['geometry']

As we are going to use the RioXarray library to mask out the data, we need to add some more information to the data. The RioXarray library is a library that allows to manipulate geospatial data in xarray datasets. Underneath it uses the rasterio library that is a library built on top of GDAL.

We need first to specify the coordinate reference system and the nodata value. Both information can be found in the metadata of the data but we need to reinforce it so that RioXarray can use it.

snow_cloud.rio.write_crs("EPSG:32632", inplace=True)
snow_cloud.rio.set_nodata(np.nan, inplace=True)

<xarray.DataArray 'stackstac-5ce38048846e4f377ce5994c6aa7ece1' (time: 101,
                                                                y: 1708, x: 1367)>
dask.array<where, shape=(101, 1708, 1367), dtype=float64, chunksize=(1, 1708, 1367), chunktype=numpy.ndarray>
Coordinates: (12/53)
  * time                                     (time) datetime64[ns] 2019-02-01...
    id                                       (time) <U24 'S2B_32TPT_20190201_...
    band                                     <U6 'scl'
  * x                                        (x) float64 6.538e+05 ... 6.811e+05
  * y                                        (y) float64 5.203e+06 ... 5.169e+06
    s2:degraded_msi_data_percentage          int64 0
    ...                                       ...
    raster:bands                             object {'nodata': 0, 'data_type'...
    common_name                              object None
    center_wavelength                        object None
    full_width_half_max                      object None
    epsg                                     int64 32632
    spatial_ref                              int64 0

Let’s clip the snow_cloud object using the catchment geometry in the UTM32N coordinate reference system.

snowmap_clipped = snow_cloud.rio.clip([geom_utm32])

It’s time to persist the data in memory. We will use the persist method to load the data in memory and keep it there until the end of the analysis.

clipped_date = snowmap_clipped.persist()

Visualize data#

We will use the hvplot library to visualize the data. The library allows to visualize data in xarray datasets. It is based on the holoviews library, which is a library that allows to visualize multidimensional data. To visualize the data on a map, we need to specify the coordinate reference system of the data. The data is in the UTM32N coordinate reference system (EPSG:32632). This will allow the library to project the data on a map. More info on the hvplot library can be found here: https://hvplot.holoviz.org/

clipped_date.hvplot.image(
    x='x',
    y='y',
    groupby='date',
    crs=pyproj.CRS.from_epsg(32632),
    cmap='Pastel2',
    clim=(-1, 2),
    frame_width=500,
    frame_height=500,
    title='Snowmap',
    geo=True, tiles='OSM')

Compute statistics#

Based on the original notebook, see the Calculate Catchment Statistics section.

We are looking at a region over time. We need to make sure that the information content meets our expected quality. Therefore, we calculate the cloud percentage for the catchment for each timestep. We use this information to filter the timeseries. All timesteps that have a cloud coverage of over 25% will be discarded.

Ultimately we are interested in the snow covered area (SCA) within the catchment. We count all snow covered pixels within the catchment for each time step. Multiplied by the pixel size that would be the snow covered area. Divided the pixel count by the total number of pixels in the catchment is the percentage of pixels covered with snow. We will use this number.

Let’s get the number of pixels of total, clouds, snow across the catchment.

# number of cloud pixels
cloud = xr.where(clipped_date == 2, 1, np.nan).count(dim=['x', 'y']).persist()

# number of all pixels per each single date
aot_total = clipped_date.count(dim=['x', 'y']).persist()

# Cloud fraction per each single date expressed in % 
cloud_fraction = (cloud / aot_total * 100).persist()

# Visualize cloud fraction
cloud_fraction.hvplot.line(title='Cloud cover %', ylabel="&") * hv.HLine(25).opts(
    color='red',
    line_dash='dashed',
    line_width=2.0,
)

We are going to get the same information for the snow cover.

snow = xr.where(clipped_date == 1, 1, np.nan).count(dim=['x', 'y']).persist()

snow_fraction = (snow / aot_total * 100).persist()

# visualize snow fraction
snow_fraction.hvplot.line(title='Snow cover area (%)', ylabel="%")

# mask out cloud fraction > 30% 
masked_cloud_fraction = cloud_fraction < 30

snow_selected = snow_fraction.sel(date=masked_cloud_fraction)

snow_selected.name = 'SCA'
snow_selected.hvplot.line(title="Snow fraction")

Let’s compare the date with the discharge data.

discharge_ds = pd.read_csv('../data/ADO_DSC_ITH1_0025.csv', sep=',', index_col='Time', parse_dates=True)

Let’s refine a little bit the data so that we can compare it with the snow cover data.

start_date = pd.to_datetime("2019/02/01")
end_date = pd.to_datetime("2019/06/30")
# filter discharge data to start and end dates
discharge_ds = discharge_ds.loc[start_date:end_date]

discharge_ds.discharge_m3_s.hvplot(title='Discharge volume', ylabel='Discharge (m$^3$/s)') * snow_selected.hvplot(ylabel='Snow cover area (%)')  

Conclusion#

In this analysis, we have comprehensively examined the features, capabilities, and limitations of two prominent geospatial data processing frameworks: OpenEO and Pangeo. OpenEO offers a unified API that simplifies the process of accessing and processing earth observation data across various backends, allowing users to interact with different data sources seamlessly. Its standardized interface is a strong asset, making it accessible to a wide range of users, from researchers to application developers.

On the other hand, Pangeo excels in facilitating big data geoscience. Its robust ecosystem, built around existing Python libraries like Dask and Xarray, makes it a powerful tool for large-scale data analysis and visualization. Pangeo’s community-driven approach and open-source nature foster collaboration and innovation, promoting a dynamic and adaptable framework.

Each platform has its own set of advantages and constraints. OpenEO simplifies interoperability and enhances accessibility, making it particularly beneficial for users who wish to work with diverse data sources without delving deeply into the complexities of each backend. Pangeo, with its emphasis on leveraging existing Python tools and libraries, is particularly potent for those comfortable with Python and seeking to perform extensive, scalable analyses.

Choosing between OpenEO and Pangeo ultimately depends on the specific requirements and constraints of a project. Considerations such as the user’s familiarity with Python, the necessity for interoperability across various data backends, and the scale of data processing required should guide the decision-making process.

How to exploit data on Pangeo

Contents

How to exploit data on Pangeo#

Snow and Cloud Cover Example#

Authors & Contributors#

Authors#

Contributors#

Introduction#

Relevant resources#

Import libraries#

Connect to the Dask cluster#

Client

Cluster Info

LocalCluster

Scheduler Info

Scheduler

Workers

Worker: 0

Worker: 1

Worker: 2

Worker: 3

Load data#

Catchment outline#

Satellite collections#

Search for satellite data using STAC#

Get bands information#

Load data#

Calculate snow cover#

Process snow cover data#

Mask data#

Aggregate data#

Visualize data#

Compute statistics#

Conclusion#